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SECTION  1 


INTRODUCTION 

Current  interests  in  aeroelastic  problems  in  the  transonic  range  demand 
viable  methods  for  unsteady  flow  calculations.  At  present,  there  are  two 
classes  of  numerical  methods  in  progress:  first,  computation  based  on  flow 
field  'meshes'  or  'elements',  and  second,  compulation  based  on  wing  planform 
elements.  The  method  in  the  first  class,  known  as  the  'computational  method', 
relies  mainly  on  the  finite  difference  scheme  (References  1  through  6)  and  on 
the  finite  element  scheme  (Reference  7).  Although  these  methods  in  many 
cases  provide  accurate  details  of  the  transonic  flow  field,  their  applica¬ 
tions  to  the  three-dimensional  flutter  calculation  appear  to  be  still  too 
costly  in  practice.  Besides,  there  are  some  frequency  restrictions  in  these 
methods  due  to  either  equation  simplification  (e.g.,  Reference  5)  or  require¬ 
ments  in  convergence  (e.g.,  Reference  l). 

The  methods  in  the  second  class  are  based  on  the  lifting-surface  theory 
formulaiion.  In  terms  of  the  paneling  scheme,  these  methods  require  consid¬ 
erably  fewer  elements  and  hence  arc  potentially  more  cost-effective. 

Ever  since  Kussner  (Reference 8V  formulated  the  classical  lifting  surface 
theory,  there  has  been  a  long-standing  research  effort  devoted  to  the  develop¬ 
ment  of  the  kernel  function  methods  (e.g..  Reference  0).  In  the  last  decade, 
this  method  was  advanced  and  adopted  by  various  aircraft  industries  as  a 
common  practice  for  flutter  calculation  purposes  (Reference  10).  As  the 
method  is  based  on  a  pressure  formulation,  it  has  the  advantage  that  the  cal¬ 
culation  domain  can  safely  exclude  the  upstream  and  the  wake  regimes  without 
any  loss  of  generality  (see  Figure  1).  Moreover,  the  method  is  proven  to  be 
llexlble  in  incorporating  various  modes  into  the  flutter  calculation  scheme  in 
i  fairly  wide  frequency  range.  The  scope  of  this  method,  however,  has  yet  to 
bo  explored  further  in  two  important  areas.  First,  the  method  was  largely 
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Figure  1.  Oscillatory  Transonic  Airfoil  as  a  Boundary  Value  Problem 


restricted  to  the  linearized  subsonic  and  supersonic  flow  regimes,  whereas 
its  application  in  the  transonic  flow  regime  only  until  recently  was  uncer¬ 
tain.  Second,  the  applications  of  the  kernel  function  formulation  have  been 
successfully  combined  with  the  use  of  the  pressure-mode  method  (References  11 
and  1 2) ,  whi ch ,  however,  lacks  the  flexible  and  versatile  feature  of  that  of 
the  panel  methods;  e.g.,  the  doublet  lattice  (DL)  scheme  (Reference  13), 

Although  the  pressure-mode  method  has  been  worked  out  in  the  purely 
subsonic  and  supersonic  regimes,  the  panel  method  has  been  limited  to  the 
subsonic  regime  and  its  supersonic  unsteady  counterpart  has  never  been  com¬ 
pletely  worked  out.  In  recent  years,  the  kernel  function  method  was  extended 
to  the  transonic  range  by  Tijdeman  and  Zwaan  (Reference  14),  and  by 
A.  Cunningham  (Reference  15),  These  methods  are,  however,  based  on  the  mod¬ 
ification  of  the  classical  purely  subsonic  and/or  the  purely  supersonic  ker¬ 
nel  functions  in  which  the  steady  flow  nonuniformity  was  not  formally  con¬ 
sidered.  Some  good  results  were  produced  by  these  methods  as  a  result  of  the 
partial  Inclusion  of  the  receding-wave  effects,  but  the  formal  expressions  of 
the  transonic  kernel  functions  and  the  calculation  procedures  remain  to  be 
establ ished. 

More  recently,  Liu  and  his  associates  (References  1  <>  through  19)  have 
proposed  a  Mixed  Kernel  Function  (MKF)  method  for  unsteady  transonic  flow 
calculations.  In  these  publications,  mixed  procedures  of  various  types  of 
kernel  functions  in  combination  with  the  panel  method  were  proposed.  The 
method  was  aimed  at  a  unified  three-dimensional  approach,  which  would  bridge 
the  classical  subsonic  and  supersonic  lifting  surface  methods  through  the 
transonic  regime,  including  the  effect  of  the  oscillatory  shock  wave.  The 
only  provision  of  the  MKF  approach  is  that  the  mean  flow  structure  in  the 
proximity  of  the  planform  must  be  supplied  by  either  experiments  or  other 
means  of  computations. 

To  tackle  the  three-dimensional  problem  as  such,  however,  requires  a 
firm  understanding  of  the  basic  MKF  procedure.  Hence,  it  is  the  purpose  of 
this  report  to  Investigate  the  basic  Issues  related  to  the  MKF  procedure  for 
two-dimensional  flow.  The  present  work  consists  of  two  goals.  On  one  hand, 
the  transonic  acceleration  potential  equation  is  derived  and  the  transonic 
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kernel  functions  in  various  flow  regions  are  established.  The  formulation 
of  the  MKF  method  is  then  followed  by  the  numerical  evaluations  of  these 
transonic  kernel  functions  and  the  application  of  the  local  linearization 
procedure  to  these  functions.  On  the  other  hand,  thorough  investigations  of 
various  discrete-element  types  of  paneling  schemes  are  required.  These 
amount  to  selecting  a  promising  candidate  from  the  DL  scheme,  the  Constant 
Pressure  Panel  (CPP)  scheme,  and  the  Linear  Pressure  Panel  (LPP)  scheme,  for 
a  unified  paneling  scheme  applicable  throughout  the  regimes  of  subsonic, 
transonic,  and  supersonic  flows.  It  is  considered  a  challenging  task,  as  no 
such  unified  paneling  scheme  has  been  found  workable  for  unsteady  flow  cal¬ 
culations  in  the  past.  In  later  sections,  it  is  demonstrated  that  our  LPP 
scheme  is  indeed  a  unified  method  suitable  for  all  linearized  flow  regimes. 

The  success  of  this  scheme  really  depends  on  the  correct  choice  of  the  con¬ 
trol  point  location.  Then  the  local  linearization  procedure  is  combined 
with  the  LPP  scheme  to  account  for  the  nonuniformity  of  the  steady  mean  flow 
(it  is  called  'nonlinear'  flow  hereafter;  specifically,  we  refer  to  it  as 
transonic/subsonic  flow,  transonic/supersonic  flow,  and  nonlinear  sonic  flow). 

To  demonstrate  the  transonic  LPP  method,  calculation  examples  are  pre¬ 
sented  for  pitching  and  flapping  NACA64A006  airfoils  in  subsonic  flow,  for 
the  pitching  Guderley  airfoil  and  parabolic-arc  airfoils  in  sonic  flow,  and 
for  the  pitching  wedges  and  parabolic-arc  airfoils  in  supersonic  flow.  Both 
purely  linearized  flow  results  and  nonlinear  results  were  compared  with  those 
calculated  by  various  methods.  A  MKF  procedure  was  performed  on  a  Guderley 
airfoil  for  the  case  of  shockless  smoothed  sonic  flow.  Also,  based  on  the 
proposed  MKF  concept,  the  oscillatory  shock  patching  scheme  was  derived  for 
the  case  of  Eckhaus-Landahl ' s  aileron  buzz  model.  Finally,  assessments  of 
and  recommendations  for  the  present  method  are  given. 


A 


SECTION  2 


BASIC  FORMULATION 


The  total  time-dependent  velocity  potential  0(x,y,z,t)  can  be  written  as 


■  U^x  +  t(x,y,z,t)  (2.1) 

where  the  barred  symbols  indicate  the  true  physical  quantities  and  being 
the  freestream  velocity.  The  perturbed  potential  i  is  normalized  by  (U^c) 
and  is  split  into  a  steady  and  unsteady  component,  i.e.. 


l/iUjc)  =  *(x,y,z,t)  =  *(x,y,z)  +  cp(x,y,z,t)  (2.2) 

where  (x,y,z,t)  =  (x/c,  y/c,  z/c ,  tU—/c)andc  is  the  physical  chord  length. 
For  transonic  flow,  the  first  term,  $,  is  governed  by  the  steady- flow  non¬ 
linear  small  disturbance  equation  (e.g.,  Reference  20)  and  the  second  term, 

9,  is  governed  by  the  oscillating  flow  equation  (e.g.,  Reference  21). 


where 


09  =  V  (#  9  ) 

®  XX  X 


(2.3) 


□  =  £  +  X _ 

dy2  3z2 

A  =  (l  -  M  2)  —  -  2ikM  3  +  ksM  3 

»  '  a  ®  ax  ® 

ax 

A  “^*tt 

9  *  9e 


(2.4) 
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ai><*  ^x  S  lx  an<*  V»  s  W  +  1)MW8.  The  parameters  M^,  y,  and  k  are  the 
freestream  Mach  number,  the  specific  heat  ratio  of  the  gas,  and  the  reduced 
frequency  (defined  as  k  ■  U)c/(2U— ),  where  c  is  the  chord  length,  and  u>  is  the 
circular  frequency  of  oscillation). 


2. 1  Transonic  Acceleration  Potential 

The  acceleration  potential,  Y,  is  related  to  the  fluid  pressure  and 
density  as 


Y  - 


(2.5) 


where  F(t)  is  some  function  of  time.  According  to  Bernoulli's  equation,  Y 
can  be  written  in  terms  of  the  total  potential,  fl(=  O/XU^c)),  i.e., 

Y  .  Q  +  (VO  •  vn  -  l)/2  (2.6) 

Substituting  Equations  (2.1)  and  (2.2)  into  Equations  (2.5)  and  (2.6)  and 
keeping  terms  of  comparable  order  to  those  in  Equation  (2.3),  we  obtain 

¥  “  *x  +  $t  +  (1  +  ♦x)^  (2,7) 


The  acceleration  potential  can  be  split  into  two  components,  i.e., 


Y 


where 


L  ♦  -  A<P 
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A 


o 


3_ 

3x 


A  -  U  +  O  +  ♦„)  |j 


(2.8) 
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It  is  seen  that  the  transonic  acceleration  potential,  contains  an  additional 
term,  ♦x,to  the  one  commonly  used  by  the  purely  subsonic  and  supersonic  flow 
analysis. 

2.2  The  Transonic  Acceleration  Potential  Equations 

Applying  operator  A  to  Equation  (2.3),  making  use  of  Equation  (2.8),  and 
keeping  only  the  lowest  order  terms  result  in 

-  V.£(*x*x)x  +  (^xx‘Px)x  +  “^xU  +  *x  •  °*x  (2‘9) 

Furthermore,  if  we  assume  that  ♦  is  a  one-dimensional  function  depending 
only  on  x,  the  last  term  in  Equation  (2.9)  becomes  (p^A^.  After  some  rear¬ 
ranging,  Equation  (2.9)  can  be  further  simplified  to  the  following  form,  i.e., 

B  a*  +  *  +  *  +  A*  +  B*  «  -ikcp  •  F  (2.10) 

o  Txx  Tyy  Tzz  x  T 


where 

6  3  =  1  -  M  3  -  X 

o  ® 

A  »  -2r  -  ik(2M  3  +  X) 

CD 

*< 

B  -  -  (y  +  2)Ao]  -  po  +  2ikro  (2.10 

F  -  u  +  ka(v  +  2)X  -  2ikr 

.  o  o  o 

Parameters  X,  X  0,  T,  T0 ,  n,  and  n0  are  all  steady  mean-flow  coefficients 
defined  as 

'  X  -  vm4>x  -  (Y  +  DX0 

.  r  -  .  (y  +  nro 

»*  •  Vxxx  “  V«M0/(VH»8  +  D  (2.12) 


J- 
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It  is  seen  that  Equation  (2.10)  contains  two  dependant  variables,  f 
and  <p,  as  related  by  Equation  (2.8).  Hence,  Equation  (2.10)  Is  In  essence 
an  lntegro-dlfferentlal  equation  in  terms  of  Jr  only.  To  simplify  the  equa¬ 
tion,  we  assume  the  mean  flow  is  accelerated  at  a  fairly  slow  rate  such  that 
T0  0(k)  and  u~0(ka),  the  right-hand  side  of  Equation  (2.11)  can  thus  be 
Ignored  up  to  the  third  order  in  k,  i.e., 

8  8*  +♦  +  A*  +  B*  -  0  (2.13) 

In  the  Addendum  to  this  report,  Landahl  derived  a  transonic  acceleration 
potential  equation  (Equation  (L.  1))  with  a  more  general  approach  in  which  he 
also  recognized  the  acceleration  potential  obtained  in  Equation  (2.8).  Start¬ 
ing  out  with  the  full  potential  equation,  he  arrived  at  a  slightly  different 
low-frequency  oscillatory  equation  from  ours  (Equations  (2.13)  and  (2.11)). 

It  turns  out  that  the  difference  lies  only  in  the  coefficients  A  and  B; 
Landahl* s  expressions  read 

* 

A  -  -2r  -  ik[2Mto8  -  (~-T)\]  +  28/ (Y  -  i)ro 

< 

B  -  kaM/  -  »  -  2ikVr0  (2.14) 

When  Equation  (2.14)  is  eompared  with  Equation  (2.11),  it  is  saen  that 
the  difference  is  insignificant  between  the  coefficients  'A',  but  there  is 
some  disagreement  between  Landahl 's  and  our  coefficient  'B'.  The  latter  dif¬ 
ference  could  be  a  result  of  our  approximation  in  starting  with  the  small  dis¬ 
turbance  equation.  However,  further  investigations  on  this  issue  are  needed. 

Next,  let  us  examine  the  coefficients,  A  and  B,  represented  by  other 
theories.  Previous  work  in  References  14,  15,  and  21  all  suggested  replacing 
the  fraestream  Mach  number  by  a  local  Mach  number  in  the  DL  method,  the 
pressure-mode  method  or  the  sonic  box  method.  We  caution  that  such  a  replace¬ 
ment  only  partially  recovers  the  mean-flow  nonuniformity,  and  that  some  terms 
were  left  out.  If  one  consistently  ignores  the  0(4/)  term  in  the  approxi¬ 
mations,  the  expressions  in  the  previous  work  then  differ  from  Equation  (2.11) 


by  a  complex  coefficient  2r  -  ikX  in  A  and  a  complex  coefficient  ka(2X  +  \,) 

+  o0  -  2ikr0  in  B.  This  implies  that  a  part  of  the  mean-flow  damping  mecha¬ 
nism  was  ignored.  In  our  earlier  work  (Reference  16)  and  that  of  Dowell  (Ref¬ 
erence  23),  both  proposed  to  employ  the  velocity  potential  equation  (Equa¬ 
tion  (2.3))  as  the  governing  equation  for  f.  In  this  case,  the  differences 
from  Equations  (2.11)  are  a  complex  coefficient,  T  +  ikX  in  A  and  a  complex 
coefficient  ka(X  +  X©)  +  Mo  -  2ikr  in  B.  With  reference  to  Equation  (8a)  of 
our  AGARD  paper  (Reference  17),  it  is  seen  that  the  coefficient  A  remains 
unchanged  but  the  coefficient  B  differs  by  a  complex  coefficient,  k*Xo  +  Mo  * 
M+2ikr.  Clearly,  the  additional  terms  are  generated  as  a  result  of  intro¬ 
ducing  the  transonic  acceleration  potential. 

In  the  subsequent  calculation  scheme,  we  have  set  up  a  completely  general 
input  for  the  coefficients  A  and  B,  so  that  one  can  easily  assess  various 
degrees  of  approximation  by  judging  from  the  numerical  results  as  well. 

Other  than  Equation  (2.11),  we  often  used,  in  the  latter  calculation,  the 
following  two  sets  of  coefficients  for  expedient  calculations.  The  first  set 
is  the  AGARD  paper  formulation  in  which  only  purely  linearized  acceleration 
potential  operator  was  used,  resulting 

A  =  -2(T  +  ikM  3 )  -  ikX 
00 

B  =  kaM  2  (2.15) 

00 

The  second  set  is  the  further  simplification  of  the  above,  i.e., 

A  =  -2(r  +  IkM3) 

00 

B  ■  kaM  a  (2.16) 

CD 

Clearly,  when  T  is  a  constant,  Equation  (2.16)  is  reduced  to  Oswatitsch  and 
Maeder's  parabolic  approximation  (References  24  and  25). 
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SECTION  3 


DEVELOPMENT  OF  TRANSONIC  KERNEL  FUNCTION 

To  derive  the  kernel  function,  it  is  necessary  to  integrate  Equation  (2.8) 
along  the  mean  flow  direction.  Hence,  the  integrated  transonic  velocity 
potential  in  the  subsonic  freestream  can  be  expressed  in  terms  of  |  as 

<P  .  .'MWf  _JL_  el«(*o)dXo  (3.1) 

-80 


where 


f(x) 


dx 

1  +  4> 

x 


and  X  =  x  -  § 


Establishing  the  kernel  function  by  the  above  equation  results  in  a 
rather  complicated  formulation.  However,  if  we  restrict  the  analysis  in  the 
low  frequency  range,  namely  0(k3,  kBi)>x),  Equation  (3.1)  can  be  approximated 
to  yield  the  commonly-used  kernel  function,  i.e., 


Km(X,Y,0) 


lim  ^7  /-oo  e1"X°1'm<xo.Y,Z)dx0 
z-»o  3z 


-ikX 

e 


(3.2) 


where  X  =  x-  5»Y  =  y~Tl*  and  subscripts  m  =  o,  1  and  2  representing  the 
sonic,  supersonic,  and  subsonic  cases,  respectively. 


3. 1  Basic  Solutions 

Basic  solutions,  ^m,  should  satisfy  Equation  (2.13).  There  are  several 
ways  to  derive  approximate  solutions  from  Equation  (2.13),  depending  on  the 


stage  at  which  we  apply  the  primary  step  of  the  local  linearization  procedure 
(Reference  26),  that  is,  to  assume  the  parameters  T,  and  li  to  be  constants. 
The  simplest  form  can  be  obtained  by  assuming  constant  coefficients  of  Equa- 
tion  (2.13)  as  the  start.  The  three-dimensional  basic  solutions  have  been 
derived  previously  in  Reference  18  (see  Equations  (10a),  (10b),  and  (10c)); 
they  read 


for  sonic  flow;  and 

=  Cj  e'AlXcosh(/AsR)/R 
for  transonic/ supersonic  flow;  and 

fe  =  C2e"AlX  + 


,  X  *  0, 


.  X  *  1  0  |r. 


for  all  X, 


for  transonic/subsonic  flow, 

where 

Aq  =  B/A 
Aj  =  A/2fc8 

•  Ag  »  -(A2  -  4etoaB)i/8/2fes 
R  *  (X2  +  ebara)1/8 
r8  =  Vs  +  z2 


(3.3) 


(3.4) 


(3.5) 


(3.6) 
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and  Cq ,  ,  and  are  the  associated  constants.  For  the  case  of 

two-dimensional  flow,  which  is  the  case  considered,  the  basic  solutions  can 
be  obtained  from  solving  the  two-dimensional  version  of  Equation  (2.13)  or 
taking  the  integral  transform  of  Equations  (3.3)  through  (3.5)  in  the  spanwise 
direction,  i.e.,  the  y-direction.  These  solutions  read 


^o  -[A°x  +  ^r] 

=  /X  * 


X  a  0, 


(3.7) 


for  sonic  flow 


♦x  =  CxeiabxJ0(oR)  ,  X  a  Bz  (3.8) 

=  0,  otherwise, 

for  transonic/supersonic  flow 

fe  =  c8elQ‘bXH,(a,(oR)  ,  for  all  X,  (3.9) 

for  transonic/subsonic  flow. 

where  J0  and  Hq  are  the  zeroth  order  Bessel  function  and  the  zeroth-order 
Hankel's  function  of  the  second  kind,  respectively;  both  functions  contain 
complex  arguments,  namely  the  complex  parameters,  a. 

The  complex  arguments  read 


„  .  iA  1A 

2  ft)*  2xaoa 

"te*  -  -ft,3  (3.10) 
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for  transonic/subsonic  flow;  the  constants  co ,  cx ,  and  09  read 


(3.13) 


It  should  be  noted  that  the  sonic  solution  Equation  (3.7)  can  also  be 
derived  as  a  limiting  case  of  the  transonic/ super sonic  basic  solution,  Equa¬ 
tion  (3.8),  and  that  of  the  transonic/ subsonic  basic  solution,  Equation  (3.9). 
When  we  let  @^2or  n\,2  approach  zero,  and  meanwhile  keeping  X  nonvanishing  in 
the  latter  equations,  it  can  be  shown  that  both  expressions  lead  to  the  sonic 
solution  by  means  of  asymptotic  expansions  of  the  Bessel  function  and  the 
Hankel  function  for  the  large  argument  oR.  Also,  we  note  that  these  basic 
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solutions  (Equations  (3.7)  through  (3.9))  contain  the  purely  sonic/supersonic/ 
subsonic  basic  solutions  (i.e.,  the  linear  solutions)  as  special  cases. 

Indeed,  when  the  transonic  parameters  T,  X,  and  tt  are  put  to  zero,  the  sonic 
solution  (Equation  (3.7))  transonic/supersonic  solution  (Equation  (3.8))  and 
the  transonic/subsonic  solution  (Equation  (3.9))  reduce  to  the  purely  linear¬ 
ized  solutions  of  Rott  (Reference  27),  of  Miles  (Reference  28),  and  of 
Posslo  (Reference  29),  as  expected.  Hence,  it  is  appropriate  to  describe 
these  basic  transonic  solutions  as  the  ' transonic/supersonic '  solution  or 
the  'transonic/subsonic'  solution,  thus,  implying  that  they  are  general  enough 
to  cover  the  transonic  flow  range  as  well  as  the  linearized  flow  range. 

3. 2  Transonic  Kernel  Functions 

The  derivation  of  the  transonic  kernel  funct:ons,  based  on  the  previous 
basic  solutions,  is  an  essential  step  towards  the  MKF  formulation.  During 
the  course  of  our  study,  the  transonic/subsonic  kernel  function  was  first 
derived,  followed  by  the  transonic/supersonic  kernel  function,  and  lastly  by 
the  sonic  kernel  function.  In  the  subsequent  development,  we  will  present 
these  kernel  functions  In  this  study  sequence. 

•  Transonic/Subsonic  Kernel  Function  (m  =  2) 

Substituting  Equation  (3.9)  into  Equation  (3.2),  and  making  use  of  the 
governing  equation  (2.13)  results  in 

-  i  kX 

Kb(X,z)  -  lim  Ia  •  e  (3.14) 

z-*o 


where 


/  2L  H,,s,<oR>«lkx°d*0 

-00  Bz 


80  el(Qb  +  k)XoHj(a1(aR)  +  igba(a0  +  k)l^<a'(aR) 

+  (V(«o  +  k)a  -  x2]  /  ei(a°  +  k)x°H)<3'  (aR)dxo 
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(  3 ) 

and  Hj  is  the  first-order  Hankel  function  of  the  second  kind.  The  last 
integral  can  be  reduced  to  a  simpler  integral  by  means  of  Laplace  transform. 

As  the  first  two  terms  vanish  in  the  far  field,  the  Iq  expression  is  further 
simplified.  Consequently,  after  the  limiting  process,  we  obtain  the  transonic/ 
subsonic  kernel  functions,  i.e.. 


MX.O)  -  ca 


eb3crx 


jYy  Hx(a)(a|  X|)  +  ieb2<ab  +  k)n>(ai<a|  X  | ) 


i(a0  +  k)X 


+  [fea(ob  +  k)a  -  x8] 


[(•  ifeMi) 


-  ikX 


(3.15) 


where 


a*  =  ia 


s  =  i  v/(Qo  +  k)3  -  a2 


l  S  =  i(afc  +  k)  +  s  (3.16) 

Several  special  cases  can  be  reduced  from  Equation  (3.15).  In  the  case 
of  parabolic  approximation  (Reference  2U)  where  X  =  p,  =  0,  the  steady  kernel 
function  can  be  obtained  by  setting  k  approach  zero.  Thus,  the  kernal  func¬ 
tions  read 


X* (X,0) 


Sgn(X)  *  Kx |  X  hK0^~  |  X 


(3.17) 
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where  Kx  andK^  are  modified  Bessel  functions  of  the  first  and  the  zer.oth 
order,  respectively.  Furthermore,  letting  f  approach  zero  results 


e. 

MX.O)  -  2X 


(3.18) 


This  is  exactly  the  purely-subsonic  steady  kernel  function. 

In  the  case  of  unsteady  purely-subsonic  kernel  function,  we  take  the 
limit  in  T  approaching  zero  for  Equation  (3.14).  The  result  reads 


K§  (X,0)  -  -  jp 


ikxo 


e  *  •  [  iMtiSgn(X)H1<a)(K  |  X  |)  -  |  X  |) 


21P-  +  rX 

—  H— ) + ‘V* 


( •) 


lkxp 

c 

(  K  I  Xo  I  )e  dxo 


-ikX 


(3.19) 


where  <  -  X/6  ,  and  Sgn(X)  -  1  if  X  >  0  and  Sgn(X)  -  -1  if  X  <  0. 


This  is  precisely  the  expression  of  Possio's  integral  equation,  which  was  derived 
by  Watkins,  Runyan,  and  Woolston  starting  from  a  three-dimensional  kernel 
function  (Reference  30,  Equation  (B. 18)). 

•  Transonic/Supersonic  Kernel  Function  (m  -  1) 

Based  on  the  nature  of  the  supersonic  wave  propagation,  the  domain  of 
Interest  for  the  transonic/supersonic  kernel  function  is  confined  to  (mz,X) 
rather  than  (-»,X).  Hence,  the  kernel  function  reads 

MX,0)  -  11m  Ij  •  e‘lkX  (3.20) 

z-*o 
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where 


I «  ■  JX  elk^  . 


dza  mz 


tld*o 


/  c,.1'*  *  k^X°Jo(oR)H(x0-  n,2>d„, 
dz  mz 


(3.21) 


Applying  the  differentiation,  making  use  of  Equation  (3.2),  and  finally  taking 
the  limiting  process  yields 


MX,0)  -  cl(noS 


Hk+d> 


ffi(X)  +  H(X)J0(OX)  - 


oH(x)J1 (aX) 


"b 


+~v 

\"b3  m o 7 


le  ikXf  H(*o>e  J0(axo)dxo 


•ikpx 

mo? 


(3.22) 


where 


ko/ny,2  =  -  (at>  +  k) 


6(X)  is  Dirac's  delta  function 


H(X)  is  the  Unit  Step  function. 


Various  special  cases  of  Equation  (3.22)  will  be  discussed  here.  For 
parabolic  approximation,  X  =*  b  »  0,  the  steady  kernel  function  in  this  case 
reads 


MX.O)  -  cim_ 


fi(X) 


+  JL  H(X)  •  x\ 

V  / 


(3.23) 


where  m  2  -  M  3  -  1 . 
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As  T  approaches  zero,  the  above  becomes 

^(X.O)  -  c^mj  •  6(X) 

On  the  other  hand,  for  purely  supersonic  flow  case,  where  T  =  \  =  u  =  0,  the 
unsteady  kernel  function  Equation  (3.22)  reduces  to 


f  -ikM.3 


MX,0)  -  Clm 


m_ 


6(X)  +  H(X)J0(«cX) 

m  2 


H!®  H(X)Jj  (*cX) 


m_ 


1.3  i  lv 

+ -  e"1  J  H(xq  )e  JoCicxoJdXo 


•ikxp 


m_ 


.  kM_  ,  2tt 

where  k  =  — 2  and  c,  =  — 
a  1  m 


(3.24) 


The  above  equation  checks  identically  with  Watkins  and  Berman's  supersonic 
kernel  function  (Equation  (49),  Reference  31). 


•  Sonic  Kernel  Function  (m  ■  0) 

The  Sonic  Kernel  function  derived  here  is  meant  to  be  the  transonic 
kernel  function,  for  =  0  (or  M  *•  1.0),  thus  including  the  transonic  thick¬ 
ness  effect.  Inserting  the  basic  solution  Equation  (3.9)  into  Equation  (3.2) 
and  making  use  of  the  governing  sonic  equation,  Equation  (2.13)  yields 


Ko(X,0) 


CoA 


<*>X 

7T-2^ 


/X  ~  ,a 
sQtoX  dX 


■ikX 


(3.25) 


19 


where 


The  sonic  kernel  function  can  be  derived  from  either  Equation  (3.15)  or 
Equation  (3.22).  It  can  be  shown  that  Equation  (3.24),  in  fact,  is  the  sonic 
limiting  case  of  these  subsonic  and  supersonic  kernel  functions.  The  limiting 
procedures  depend  primarily  on  the  validity  of  the  asymptotic  expansion  of 
the  Hankel  function  and  the  Bessel  function,  provided  that  the  domain  of  'X1 
(of  the  argument  cX)  be  excluded  from  zero. 

For  the  case  of  purely  sonic  flow,  we  let  T  =  X  ■  u  m  0,  and  let  »  1.0, 
Equation  (3.25)  then  is  reduced  to  the  same  form  as  the  equation  itself  except 
that  the  coefficients,  c©  and  a0 ,  are  simplified  to 


, .  /rrik 
CO  V  — 

~  ik 

00  “  2  (3.26) 


The  purely  sonic  kernel  function  here  checks  identically  with  that  derived  by 
Watkins,  et  al.,  (Equation  (B.23),  Reference  30). 
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SECTION  A 


PANEL  METHODS  FOR  PRESSURE  EVALUATION 

The  objective  of  the  present  formulation  Is  to  make  use  of  the  kernel 
function  approach  for  solving  the  downwash  Integral  equation,  l.e., 

l  1 

V(x)  "  2^  J*  Km<X.0>P<5)dC  (4.1) 

where  l  Is  the  total  normalized  chord  length  (normalized  by  half-chord  length, 
l.e.,  I  •  2)  and  the  pressure  function  Is  the  solution  sought, defined  as 

p<?)  "  [plower(?)  -  pupper(?)l  /(lp«U.a)  (4‘2) 

Given  the  downwash  function  W(x),  one  Is  required  to  solve  for  p(§). 

The  solution  techniques  for  pressure  In  the  kernel  function  approach  can  be 
broadly  categorized  In  two  classes,  namely,  the  pressure-mode  method  (e.g.. 
References  11  and  12)  and  the  panel  method  or  the  discrete  element  method, 
e.g.,  References  10,  13,  and  1A.  The  general  approach  Is  to  prescribe  an 
approximate  solution  p(5)  In  the  following  form,  l.e., 

N 

P(?)  -  E  P,US>  <4. 3) 

1-1 

where  5  and  the  subscript '1'  represent  the  sending  point  and  the  sending 
panel.  In  the  pressure-mode  method,  further  knowledge  of  f^  is  needed,  a 
priori.  For  the  panel  methods,  f^  can  be  expressed  more  freely  in  various 
ways,  depending  on  the  sophistication  of  the  paneling  scheme  as  one  desires. 

There  are  several  panel  methods  now  being  widely  used.  The  simplest 
method  is  the  doublet  lattice  (DL)  method  (e.g.,  Reference  13).  It  is  a 
commonly-used  tool  for  flutter  analyses.  But  the  method  Is  restricted  to  the 
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subsonic  flow  regime  due  to  Its  Inherent  elliptical  characteristic  of  the 
collocation  scheme  adopted.  The  constant  pressure  panel  (CPP)  method  was 
first  proposed  by  Woodward  (Reference  32)  and  was  extended  to  cope  with  the 
unified  unsteady,  purely  subsonic/supersonic  three-dimensional  flow  cal¬ 
culations  by  Winther  (Reference  33).  While  the  unsteady  part  of  the  suosonic 
calculation  worked  out  successfully,  the  unsteady  part  of  the  supersonic  flow 
calculation  did  not.  This  then  led  us  to  investigate  the  applicability  of 
the  Linear  Pressure  Panel  (LPP)  method  for  a  unified  unsteady,  subsonic/ 
transonic/supersonic  flow  calculations.  The  LPP  method,  also  known  as  the 
elementary  vortex  distribution  (EVD)  method,  was  first  proposed  by  Shen,  et  al., 
(Reference  34)  for  incompressible  steady  flow  calculation  for  lifting  wings 
with  jet-augmented  flaps.  Apart  from  the  unsteady  subsonic  work  by  Fromme, 
et  al.,  (Reference  35),  no  attempt  was  made  to  unify  such  an  approach  to  the 
supersonic  flow  regime.  In  what  follows,  we  shall  discuss  the  DL,  the  CPP, 
and  the  LPP  paneling  schemes,  with  special  emphasis  on  the  development  of  the 
LPP  method  for  unified  unsteady  transonic  flow  calculations. 


4.1  POWNWASH  MATRIX  EQUATION  AND  PRESSURE  FUNCTIONS 

The  downwash  integral  equation  (4.1)  in  combination  with  Equation  (4.3) 
can  be  written  in  terms  of  the  panels  scheme  as 


W(x) 


2tt 


N 

E 

i-1 


k(x  -  5)f±<5)d5 


(4.4) 


where 

N  *  total  number  of  panels 

Pj  »  pressure  at  the  sending  point  !^(DL),  or  pressure  magnitude  on  the 
panel  (CPP),  or  pressure  at  the  leading  edge  of  the  ith  panel  (LPP) 
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Doublet  Lattice 
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*  ♦ 

x 


f.. 

i 


***  ft** 


■l+l 


Constant  Pressure  Panel 


Linear  Pressure  Panel 


Figure  2.  Sending  and  Receiving  Points 

(or  Panels)  for  DL,  CPP,  and  LPP  Methods 


Now,  let  W  =  W(xj )  where  is  the  location  of  the  receiving  point  (or  called 
the  control-point)  associated  with  the  jth  panel.  In  terms  of  the  matrix  element 
form.  Equation  (4. 4)  reads 

wj  '  D)i  -  pi  «-5> 


where 


and 


’i+l 


Dji  '  75  J.  «*j  ‘  «»,«>« 

5i-l 


i  *  l,  ....  N 

J  ■  1.  ....  N 

N  being  some  positive  integer 


(4.0) 
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In  the  ith  panel,  the  pressure  fuction  f ^ ( 5 )  is  prescribed  by  various  models 
according  to  the  panel  schemes  used. 


For  the  DL  method, 

fi<5>  *  «(?  -  ?i)  • 

where 

■  length  of  the  ith  panel 

4(5)  *  Dirac's  delta  function 
For  the  CPP  method, 

f i (5)  -  H^)  -  H(5i+1> 


(A. 7) 


(A. 8) 


where 

H(5)  *  the  unit  step  function 

For  the  LPP  method,  the  ith  regular  panel  element  has  the  pressure  function 


(?i+i  -  5)/(S1+1  -  §,)  .  5,  *  S  *  51+1 

ft<5)  -  < 

(5  -  ?i.1)/(5i  -  S^j)  ,  Slml  *  5  *  ?!  (4.9) 


p 


4. 2  The  Kernel  Matrix  Element  < 

Accordingly,  the  kernel  matrix  element  D^,  Equation  (4.5)  can  be  expressed, 
in  terms  of  Equations  (4.7)  through  (4.9),  as  follows. 

For  the  DL  method 


Dji  >  ~  Ke(xj  -  li>^i  (4.10) 

i  =  1,  ....  N 

J  “  1*  ...»  N 


For  the  CPP  method 


#i+l 


D)i  ■  is  S,  Vi  ■ ?)d5 

5i 


For  the  LPP  method 


*i+l 


D)i  *  15  /.  V*j  -  5>fi<5)d5 

5i-l 


(4.11) 


(4.12) 


Notice  that  the  domain  of  integration  for  the  ith  panel  is  C?j»  for 

the  CPP  method  and  is  for  the  method.  Clearly,  these  three 

methods  differ  from  each  other  in  their  ease  of  application. 

The  DL  method  requires  no  integration  for  obtaining  Dj^,  hence  is  the 
simplest,  whereas  the  CPP  and  the  LPP  method  require  the  integration  proce¬ 
dure.  In  fact,  in  terms  of  the  integration  scheme,  the  LPP  method  is  no  more 
than  a  'higher-moment'  integral  of  the  CPP  method;  hence  the  integration 
procedure  of  the  LPP  method  is  slightly  more  complicated.  On  the  other  hand, 
it  is  this  'higher-moment'  characteristic  of  the  LPP  method  that  provides  a 
linear  distribution  of  the  calculated  pressure,  which  is  a  definite  refinement 
of  the  CPP  method.  The  common  characteristics  of  the  CPP  and  the  LPP  methods  are 
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that  the  edge  behaviors  (including  leading  edge  singularities,  trailing  edge 
behaviors  and  hinge-point  singularities.  Reference  36)  and  the  shock-jump 
behaviors  (e.g.,  References  37  and  38)  to  a  large  extent,  can  be  incorporated 
into  their  kernel  matrix  formulation,  whereas  the  DL  method  does  not  have 
this  advantage.  However,  these  edge  behaviors  must  be  derived  from  analyti¬ 
cal  means.  Such  types  of  edge  behavior  and  shock-jump  behavior  have  not  yet 
been  formally  derived  for  airfoils  in  the  unsteady  transonic  regime  (now  with 
the  exception  of  Reference  38).  For  this  reason,  in  the  subsequent  calcula¬ 
tion,  we  use  either  the  known  purely  subsonic  leading  edge  singularity  or, 
for  the  most  part,  set  up  the  calculation  scheme  without  the  incorporation  of 
the  edge  behaviors.  In  the  following  paragraphs,  however,  pressure- function 
models  are  set  up  at  various  elements  in  different  flow  regimes.  These  pres¬ 
sure  models,  when  properly  assigned,  indeed  provide  correct  edge  behaviors. 
This  can  be  observed  from  the  computed  numerical  results  (see  Section  5  and 
figures).  However,  one  is  cautioned  that  the  'singularity-like*  leading  edge 
(for  subsonic  or  sonic  flow)  or  the  hinge  point  behaviors  obtained  are  really 
not  singular,  if  one  takes  up  the  present  element  modeling  for  f^(?)  in 
paragraph  4.3.  We  achieve  these  'singularity-like'  edge  behaviors  by  refining 
the  element  size  in  the  neighborhood  of  these  edges.  The  usual  practice  in 
our  computation  scheme  is  to  distribute  about  five  to  ten  elements  within 
five  percent  chord  length  in  the  edge  neighborhood.  In  so  doing,  one  obtains 
more  enhanced  singularity- like  behavior.  What  really  happens  in  this  outcome 
is  the  following.  Although  the  f^(5)'s  are  not  singular,  and  the  correct 
f^(S)'s  should  be  singular  (but  are  integrable  singularities),  always 
remains  regular  in  either  case.  Also,  the  value  of  would  amount  to  the 
same  order  of  magnitude  locally  for  both  cases.  This  then  in  turn  assures 
the  singularity- like  behavior  in  the  pressure  function  at  these  edges. 

In  what  follows,  we  will  restrict  our  considerations  to  the  LPP  method, 
as  it  constitutes  the  main  body  of  the  present  MKF  calculation  procedure. 
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4.3  The  Transonic  LPP  Model  of  tj(§) 

The  success  of  the  LPP  method  depends  mainly  on  the  correct  modeling  of 
the  pressure  function,  hence  fi(5).  Here,  regular  elements  and  the  irregular 
elements  will  be  described.  The  regular  elements  consist  of  the  usual  reg¬ 
ular  panel  element,  the  leading-edge  element,  and  the  hinge-point  elements. 
The  irregular  elements  consist  of  trai l ing-edge  elements  and  the  shock-point 
element.  The  last  element  will  not  be  included  in  the  following  description. 


•  Regular  Element 

The  typical  regular  element  is  shown  in  Figure  3a  and  the  pressure 
function  is  defined  as  Equation  (4. 9). 

•  Leading-Edge  Element 


Figure  3b.  Leading-Edge 


^.  =  0  (i.e.,  the  leading  edge) 


5a  -  5 

MO)  =  1  (A. 13) 


E 1 ement 


As  mentioned  earlier,  it  is  seen  that  the  present  modeling  of  fj(5) 
does  not  provide  the  leading-edge  singularity. 

•  Hinge-Point  Element 

The  pressure  function  behind  the  hinge  point  of  the  oscillating  flap 
assumes  the  same  pressure  function  as  that  of  the  leading-edge  ele¬ 
ment,  i.e.,  Equation  (4.13),  except  that  '  is  now  the  hinge-point 
location.  Ahead  of  the  hinge  point  ,  the  pressure  function  then 
assumes  a  mi rror  image  of  ft(5)  along  the  hinge  line  5  *  5i  •  It  should 
be  remarked  that  although  the  strength  of  the  edge  singularities 
varies,  depending  on  the  location  (at  the  leading  edge  or  at  the 
hinge  point),  the  local  Mach  number,  and  airfoil  geometry,  they  are 
usually  integrable  singularities,  even  for  the  cases  in  oscillatory 
flow.  Consequently,  the  kernel  matrix  here  remains  regular 
regardless  of  the  edge  behavior. 
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•  Supersonic  Trai l ing-Edge  Element 


Figure  3e.  Supersonic  Tral l lng-Edge  Element 


5N  _  2 


ll.e.,  the  trailing  edge,) 


fN<?)  ■ 


<5n+i  '  «>'<Vi  -  2) 


<5  -  Vl>/(2  -  Vl>  Vi5?82 


(4.16) 


It  Is  necessary  to  clearly  define  the  regular  subsonic,  sonic  and 
supersonic  elements  at  this  stage.  Given  the  mean  flow  pressure  (or  local 
Mach  number),  we  can  measure  the  local  value  of  0 q  at  the  mid  chord  of  the 
element.  The  panel  element  is  then  defined  as  the  subs  >nic  element,  super¬ 
sonic  element,  or  sonic  element  depending  on  whether  P  is  positive,  negative, 
or  zero,  respectively.  For  the  case  of  mixed  flow  (e.g.,  Figure  3f),  two 
extra  types  of  panel  elements  are  added. 

First,  the  panel  elements  are  tailored 
to  fit  the  sonic  point  location;  this 
results  In  two  adjacent  panel  elements, 
one  on  the  upstream  side  of  the  sonic 
point,  the  other  on  the  downstream  side. 

('bus,  these  adjacent  elements  are 
defined  as  the  subsonic/sonic  element 
and  the  sonlc/supersonlc  element  (Fig¬ 
ure  If). 


Sonic  Line 
/ 

M  <  1  /  M  >  1 

/ 

/ 

/  Sonic/supersonic 
/  element 

i  i  0  <i  ■  (■)  i  ■  .  -*■  - 

Subsonic/sonic 
el oment 


Figure  3f.  Adjacent  Elements  in 
Mixed  Flow 


4.4  The  Control  Point  x^ 

It  is  crucial  to  determine  the  proper  control  point  (or  receiving  point, 
Xj)  location  for  the  LPP  elements  in  various  local  flow  regimes  in  order  to 
obtain  the  correct  pressure  distribution.  A  classical  case  is  the  Doublet 
Lattice  method  in  which  remarkable  subsonic  pressure  distributions  result 
when  the  control  point  is  located  at  the  three-quarter  chord  of  each  ele¬ 
ment  (References  13  and  39). 

As  a  result  of  our  investigation,  we  found  that  the  control-point  loca¬ 
tion  of  each  LPP  element  is  not  entirely  independent  of  the  local  flow  con¬ 
ditions.  For  regular  elements  and  irregular  trailing-edge  elements,  the 
control-point  locations  are  determined  by  the  following 


x 


J 


'  (Sj  +  5)+1>/2 
5; 


for  subsonic  elements 
and  sonic  elements 

for  supersonic  elements 


(4.17) 


This  is  to  say  that  the  proper  location  of  the  control  point  should  be 
placed  at  the  midchord  of  the  subsonic  and  the  sonic  elements,  whereas  it 
should  be  placed  at  the  leading  edge  of  the  supersonic  element.  In  the 
mixed-flow  case,  the  adjacent  subsonic/sonic  element  assumes  the  same  control 
point  location  as  that  of  the  subsonic  element.  For  the  adjacent  sonic/ 
supersonic  element,  we  place  the  control  point  at  the  midchord. 

We  are  fairly  confident  of  the  choice  of  the  control-point  locations  for 
the  former  elements.  But  for  the  adjacent  element,  our  choice  of  the  control 
point  may  not  be  the  most  appropriate,  as  shown  in  the  sonic  flow  calculation 
(Figure  31).  Hence,  further  investigation  on  the  mixed  flow  panel  elements 
is  needed. 
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4.S  Evaluation  of  Kernel  Matrix  Element  0 j  ^ 

To  evaluate  the  Integral  (Equation  4.6) 


*1+1 


‘)l  '  75  X,  V]  '  5)f,(5)<15 

*1-1 


It  Is  necessary  to  assess  the  singular  behavior  of  the  transonic  kernel  function 
Km(X).  In  the  previous  section,  we  have  shown  that  while  the  subsonic  and  the 
sonic  kernels  are  singular,  Kp(X)  ~  1/X  and  K0(X)~  1//X,  the  supersonic  ker¬ 
nel  K1(X)~6(X1  as  | X |  approaches  zero.  Consequently,  regular  functions 
F'm  (m  "  0,1,2)  are  Introduced  for  the  purpose  of  evaluation. 

For  convenience,  let  us  first  define  the  element  Integral  operator  In 
Equation  (4.b)  as 


V-I'tsX  l-l« 

*1-1 


*1+1 


(A. 18) 


Thu  s , 


a  Subsonic  Element  (m  ■  2) 


F » <xj  -  5 )  . 

'ji  ■  h  i  -  -  i  ■  fi<5>1 


-  5)  -  <Xj  -  5)K8(xj  -  ?) 

•  Supersonic  Element  (m  -  1) 


Jl 


(A. 19) 


-  5)  *  fi<S»  +"5.8<xJ 


-  5^  ,  Xj  > 


(A. 20) 


*  XJ  4  ?i-l 
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M*j  -  V 


Kx(xj  -  5)  -  — *5)  *  xj  *  5 


.  <C 


Sonic  Element  (m  -  0) 


(F0(x. -  5) 

‘.h =  ’  '*(?> 

l  -5 


'  X1  >  5i-i 


»  X1  <  5i-1 


(A. 21) 


F0<Xj  -  5>  " 


K0(x, "  5) 


■  xj 2  5 


,  Xj  <  ? 


Now,  since  F  (X)'s  are  analytic  everywhere  over  the  domain  of  interest 

in  Xj  and  5,  the  singularities  in  Dj^  involving  cases  of  m  «  2  and  0  are 

removed.  Thus,  the  functions  Fm<X)  are  subject  to  the  parabolic  curve  fit 

scheme.  The  values  of  F  (X)  are  obtained  at  three  points  of  each  element, 

m 

namely,  the  leading  edge,  the  midchord  and  the  trailing  edge.  In  the  case  of 
sonic  and  sonic/supersonic  elements  (with  control  point  located  at  the  mid¬ 
chord),  we  simplify  the  curve  fit  scheme  by  a  linear  fit  in  letting  Fm(X) 
go  through  the  leading-edge  point  and  the  midchord. 

As  an  example,  we  shall  demonstrate  the  evaluation  scheme  for  a  regular 
subsonic  element,  using  Equation  (4.9)  and  Equation  (4.19),  First,  let  us 
write  the  curve- fit  formula,  i.e., 


F,(Xj  -  5)  -  (Xj  -  5)K,(Xj  -  5) 


0  <  Xj  «  2 


A(Xj  -  5)=  +  B(Xj  -  5)  +  C 


(4.22) 
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where  A,  B,  and  C  are  determined  by  the  parabola  fit.  Now,  according  to 
Equation  (4.9)  and  Fa  (X)  given  above,  can  be  integrated,  for  the  ith 
or  the  (i+1) th  element  (denoted  as  AD^),  to  yield  (refer  to  Figure  3a) 


AD 


Ji 


f*«f(x  -  i)  +  B6l+  2C[l  +  (x  -  1 )  in  |  T}  (4.23) 


where 


x  -  <xj  -  5i>/6i 
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*  ?i-l  s  ?  5^ 
»  ?i  4  §  «  Si+1 


(4.24) 


Notice  that  now  AD^  is  a  regular  function  as  a  result  of  the  proper 
arrangement  of  Fa(X).  For  the  evaluated  formulas  of  for  the  sonic  and 

the  supersonic  elements,  we  refer  them  to  Appendix  A. 
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SECTION  5 


RESULTS  AND  DISCUSSION 


In  this  section,  numerical  results  for  all  flow  regimes  are  presented 
and  followed  by  brief  discussions.  Figures  4  through  13  present  all  the 
transonic/subsonic  results  for  thin  airfoils  (flat  plates)  and  for  a  pitch¬ 
ing  or  flapping  NACA64A006  airfoil,  using  either  the  Transonic  Doublet  Lattice 
method  (TDL2D  code)  or  the  Transonic  Linear  Pressure  Panel  method  (TLP2D 
code).  Based  on  the  unified  formulations  in  the  previous  section,  we  shall 
demonstrate  further  that  the  LPP  method  is  equally  applicable  to  the  sonic 
and  supersonic  cases,  whereas  the  DL  method  is  not.  Figures  14  through  22 
present  transonic/supersonic  results  for  pitching  thin  airfoils  (flat  plates), 
wedges,  and  parabolic-arc  airfoils.  Sonic  flow  results  are  shown  in  Figures 
24  through  31,  for  pitching  thin  airfoils,  Guderley  airfoils,  and  parabolic- 
arc  airfoils.  In  particular,  Figure  23  illustrates  the  divergence  of  the 
purely  subsonic  and  purely  supersonic  kernel  calculations,  and  Figure  31 
presents  results  of  a  study  for  the  mixed  kernel  function  procedure. 

Notice  that  in  all  the  figures,  only  the  in-phase  and  out-of-phase 
pressure  distributions,  ACp'  and  ACp",  or  their  magnitudes  and  phase  angles, 

|  ACp  |  and  $,  are  presented.  The  pressure  coefficients  are  defined  by 
Equation  (4.2)  as 

ACp'  =  Real  part  [p],  the  in-phase  pressure  and 

ACp"  =  Imaginary  part  [p"J,  the  out-of-phase  pressure. 

All  the  pressure  coefficients  are  plotted  against  the  chordwise  coordinate,  x. 
For  convenience  of  presentation,  the  range  of  x  and  that  of  5  is  changed  to 
(0,1),  this  amounts  to  changing  the  total  chord  length  l  =  1,  instead  of  2, 
as  defined  in  Section  4.  The  percent  thickness  in  all  the  plots  is  defined 
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Figure  8.  Purely  Subsonic  and  Transonic/Subsonic  Pressure  Distributions  using  the  DL  Method 
and  the  LPP  Method  for  NACA64A006  Airfoil  at  M,  =  0.825  and  k  =  0.3  with  Pitching 
Axis  at  Midchord 
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Figure  9.  Flapping  NACA64A006  Airfoil  at 
-  0.825  and  k  ■  0.062  with 
Hinge  Point  at  Three-Quarter 
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Figure  11.  Phase  Angles  for  a  Flapping  NACA6AA006  Airfoil  at  M  =  0.8  and  k  =  0.253  with 
Hinge  Point  Located  at  Three-Quarter  Chord 
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Figure  1A.  Pitching  Thin  Airfoil  (t  -  0)  at  «  1.2  and  Reduced  Frequency 
k  -  1.0  with  Pitching  Axes  at  (a)  Leading  Edge,  (b)  Midchord, 
and  (c)  Trailing  Edge 


Figure  15.  Pitching  Thin  Airfoil  (t  ■  0)  at  ■  1.15  and 
Reduced  Frequencies  at  (a)  k  =  0.2  and  at  (b) 


k  =  0.6  with  Pitching  Axis  at  the  Leading  Edge 
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Figure  24.  Pressure  Magnitudes  of  Pitching  Thin 
Airfoils  (t  =  0)  at  Me#  *  1.0  and  at 
Reduced  Frequencies  k  =0.05  (Circles) 
and  k  =  0.5  (Triangles)  with  Pitching 
Axis  at  the  Leading  Edge 
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Reduced  Frequencies  k  =  1.0  (Circles) 
and  k  =  2.5  (Triangles)  with  Pitching 
Center  at  the  Leading  Edge 
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Figure  26.  Phase  Angles  of  Pitching  Thin  Air¬ 
foils  (t  *  0)  at  Mgs  ■  1.0,  with 
Pitching  Axis  at  the  Leading  Edge, 
using  the  Sonic  Linear  Theory  (Ref 
26,  open  symbols)  and  the  Present 
LPP  Method  (solid  symbols) 
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Figure  30.  Comparison  of  Pressure  Distributions  for  a  Six  Percent 
Parabolic-Arc  Airfoil  in  Nonlinear  Sonic  Flow  (Me*  "  1.0} 
k  *  0)  at  Steady  Angle  of  Attack 


as  two  times  that  of  the  thickness  ratio,  T,  and  T  Is  defined  as  t/c  for  an 
airfoil  and  0/2  for  a  finite  wedge  where  t  Is  one-half  of  the  maximum  thick- 

A 

ness  and  0  Is  the  apex  angle. 

5. 1  Computer  Program  Description 

The  general  inputs  to  the  TDL2D  and  the  TLP2D  programs  consist  of  the 
given  flow  parameters  and  the  paneling  parameters.  The  given  flow  parameters 
are  the  freestream  Mach  number,  M*.  the  reduced  frequency,  k,  and  the  pitch¬ 
ing  axis  location,  x^ ;  the  paneling  parameters  are  the  total  number  of  panels 
assigned  along  the  chord,  the  sending  point  location  (for  TDL2D  only),  and 
the  control  point  location.  Also,  the  distribution  of  the  panels,  fine  or 
coarse,  can  be  easily  tailored  in  these  programs.  For  nonlinear  transonic 
flow  calculations,  the  transonic  parameters  and  the  steady  mean- flow  input 
must  be  considered.  The  coefficients  A  and  B  are  built  in  with  options  of 
using  Equation  (2.11),  (2.15),  or  (2.16);  the  transonic  parameters,  X,  T,  and 
H  can  also  be  assumed  different  values,  or  put  to  zero  individually.  These 
parameters  become  functions  of  the  steady  mean-flow  pressure  (see  Equa¬ 
tion  (5.2)),  when  the  methods  are  switched  to  the  scheme  of  local  linearization. 

The  steady  flow  input  is  always  given  by  other  means,  either  by  other 
theories,  numerical  methods,  or  by  experiment.  For  example,  the  steady  flow 
inputs  for  the  NACA04A006  airfoil  in  Figures  8  through  13  are  based  on 
Tijdeman  and  Schippers'  measurement  (Reference  41),  whereas  in  the  supersonic 
cases  (Figures  19  through  22)  and  sonic  cases  (Figures  27  through  31),  the 
inputs  are  based  on  analytical  formulas.  At  any  rate,  the  steady  pressure 
coefficient  CpQ  is  expressed  in  terms  of  a  polynomial,  i.e., 

Cp  -  I  a„5n  .  0  <  5  <  1  (5.1) 

o  n*0 

where  §  is  the  'sending-point'  coordinate,  N  is  an  integer  up  to  N  =  9.  Of 

course,  an's  are  obtained  by  applying  the  curve-fit  scheme  to  the  given  data. 

The  C_  's  used  for  all  figures  are  listed  in  Appendix  C.  Thus,  the  transonic 
“ o 

parameters  read 
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W=^2cp"(§)  (5.2) 
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g 

In  the  TLP2D  code,  the  parameter  ,  rather  than  M,  is  used  as  an 

indicator  of  the  type  of  the  local  flow  characteristics,  and  hence  the  type 

of  panel  element  that  should  be  applied.  Given  C_  ,  Equation  (5.1), 

Fo 

the  calculation  scheme  is  switched  to  subsonic,  sonic,  or  supersonic  type 
according  to  wherever  0q8  is  less  than,  equal  to,  or  greater  than  zero,  mea¬ 
sured  at  the  midchord  (usually  the  control  point)  of  ea,ch  element.  Such  a 
local- flow  indicator  is  specially  designed  for  the  case  of  mixed  flow,  hence 
the  mixed  kernel  function  procedure.  More  will  be  discussed  in  the  sonic 
flow  section. 

5. 2  Transonic/Subsonic  Results 

The  in-phase  and  the  out-of-phase  pressure  distributions  (i.e.,  ACp' 
and  ACp")  for  flapping  and  pitching  NACA64A006  airfoils  in  purely  subsonic 
flow  are  presented  in  Figures  4,  5,  6,  and  7.  The  purpose  here  is  to  check 
out  the  TDL2D  and  the  TLP2D  codes  against  the  linear  theory.  It  is  seen 
that  the  present  results  are  in  good  agreement  with  the  linear  theory  and 
those  by  the  finite-difference  method  (References  3,  and  42).  It  should  be 
remarked  that  in  all  these  cases,  the  results  obtained  using  the  TDL2D  code 
are  indistinguishable  from  those  in  Figures  4  through  7  using  the  TLP2D  code. 

Figures  8  through  13  present  the  transonic/subsonic  calculations  in  the 
subcritical  regime  for  NACA64A006  airfoils  in  pitching  and  flapping  motions. 

In  these  cases,  the  steady  mean- flow  inputs,  Cp^'s,  are  given  by  TiJ deman  and 
Schippers  (Reference  41),  and  are  curve-fitted  as  polynomials  (see  Appendix  C). 
Figure  8  compares  the  linear  and  the  nonlinear  results  calculated  by  TDL2D 
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and  TLP2D  codes  for  a  pitching  airfoil.  Figures  9  through  13  show  the 
nonlinear  results  based  on  TDL2D  or  TLP2D  codes  for  flapping  airfoils.  It 
can  be  observed  that  all  theoretical  methods  produce  pressure  magnitudes 
generally  larger  than  the  experimental  data  (Reference  41).  Such  a  differ¬ 
ence  could  be  due  to  neglecting  the  boundary  layer  effects  in  the  analysis. 

For  detailed  discussion  on  this  issue,  one  is  referred  to  References  14,  19, 
and  53,  In  Figure  11,  the  present  results  in  phase  angles  are  seen  to  have 
some  deviation  from  the  finite-difference  results  (References  l  and  42)  from 
the  leading  edge  to  midchord.  In  Figures  12  and  13,  various  results  in  pres¬ 
sure  magnitudes  and  phase  angles  are  plotted  using  the  TDL2D  code,  TLP2D  code, 
and  Ehler's  finite-difference  method.  Again,  it  is  seen  that  there  is  some 
disagreement  between  the  present  results  and  the  experimental  data  behind 
the  leading  edge.  Although  the  TDL2D  result  checks  well  with  the  TLP2D 

result  in  AC  ' ,  considerable  difference  is  seen  in  the  AC  "  results.  The 
P  P 

cause  of  this  difference  is  not  clear  and  will  be  subject  to  future 

investigation. 

Some  remarks  are  in  order  with  regard  to  our  experience  in  using  the 
developed  codes.  First,  as  pointed  out  earlier,  our  results  exhibit  some 
discrepancies  in  the  transonic/subsonic  calculations,  particularly  near  the 
leading  edge.  Our  experience  shows  that  these  are  caused  mainly  by  the 
leading  edge  singularity  of  the  steady  mean-flow.  This  type  of  singularity 
actually  results  from  the  i l l- formulation  of  the  small  disturbance  theory  in 
the  leading  edge  region,  be  it  linear  or  nonlinear,  due  to  the  Incompatibility 
of  the  equation  and  the  linearized  boundary  condition.  In  other  words,  the 
singularities  may  vary  in  their  'strengths'  in  different  flow  regimes,  but  once 
they  are  adopted  as  the  steady  mean-flow  inputs  for  our  transonic  parameters 
(Equation  ( 5. 2)),  discrepancies,  or  even  numerical  instabilities  may  occur. 

(In  particular,  the  most  sensitive  parameter  of  all  is  the  acceleration 
parameter,  T,  which  is  often  very  large  in  the  neighborhood  of  the  leading 
edge. ) 
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To  attain  a  stable  computation  procedure,  one  is  required  to  keep  the 
value  of  T  constant  and  let  it  be  no  greater  than  order  of  one  in  the  first 
few  panels.  Then  T  resumes  its  local  values  at  the  panel  where  the  true 
value  of  T  is  Indeed  below  one.  Strictly  speaking,  this  procedure  is  not 
altogether  legitimate.  In  the  Addendum,  Professor  Landahl ' s  general  formula¬ 
tion  of  the  phase  correction  method  (Equation  L. 4)  may  circumvent  such  a  pro¬ 
cedure.  In  his  formulation,  the  T-term  can  be  integrated  out  altogether  as 
a  result  of  a  further  delay  step  of  the  local  linearization  application.  For 
this  reason,  it  is  felt  that  some  serious  attention  must  be  given  to  Landahl 's 
phase-correction  method. 

Next,  we  need  to  describe  the  ease  of  application  of  the  TDL2D  code  and 
the  TLP2D  code.  In  all  the  calculations  using  the  TDL2D  code,  we  only  require 
equal  spacing  panels;  normally  we  set  N  »  20,  regardless  of  pitching  or  flap¬ 
ping  motion.  In  the  cases  of  using  TLP2D  code,  more  refined  panels  are  required 
in  the  neighborhood  of  the  leading  edge  and  the  hinge  points,  the  panel  num¬ 
bers  ranges  N  ■  23  to  N  =  33.  As  the  LPP  method  is  a  higher  moment  method 
than  the  DL  method  (see  Section  4),  the  total  CPU  time  for  the  same  case 
required  using  TLP2D  code  runs  usually  more  than  twice  of  that  using  TDL2D 
code. 

5. 3  Transonic/Super  sonic  Results 

To  check  out  the  LPP  method  in  the  supersonic  regime,  calculations  using 
the  TLP2D  code  are  first  performed  for  several  purely  supersonic  examples 
(Figures  14  through  18).  Computed  results  for  pitching  flat  plates  using 
TLP2D  are  compared  with  Chadwick-Platzer ' s  result  using  the  linearized  method 
of  charac teristics  (MOC)  (Reference  43)  in  Figures  14  and  IS;  also,  computed 
results  are  compared  with  Jordon's  theoretical  result  (Reference  44)  for  in- 
phase  and  out-of-phase  pressure  distributions  in  Figure  lb.  In  these  fairly 
high  frequency  ranges  (k  ■  0.5  and  0.6),  it  is  seen  that  the  present  results 
are  in  good  agreement  with  the  other  results.  Figures  17  and  18  again  show 
the  comparisons  of  the  out-of-phase  pressure  distributions  of  the  present 
method  with  those  obtained  by  the  classical  method  of  Garrick  and  Rubinow 
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(Reference  46),  by  Brix-Platzer ' s  linearized  MOC  (Reference  47)  and  by 
Snyder's  finite  difference  method  (Reference  45).  Our  results  are  in  excel¬ 
lent  agreement  with  the  results  calculated  according  to  Garrick  and  Rubinow. 
Figures  19  through  22  show  the  transonic/supersonic  calculations  in  the  non¬ 
linear  supersonic  (or  low- supersonic)  regime  for  wedges  and  for  parabolic- 
arc  airfoils.  Using  the  linear  formula  (Equation  (4.26),  Reference  54)  for 
a  wedge  (T  =  0.05)  as  steady  mean- flow  input  (see  Appendix  C),  we  obtained 
the  out-of-phase  pressures  at  two  different  frequencies.  The  comparisons  in 
Figure  19  show  that  our  result  is  in  better  agreement  with  Carrier's  exact 
theory  than  Chadwick-Platzer' s  MOC  method.  Notice  that  the  present  calcula¬ 
tion  scheme  for  the  wedge  case  makes  no  approximations,  as  the  transonic 
parameter,  X,  is  constant  and  T  and  |J<  are  identically  zero.  Figure  20  shows 
the  nonlinear  results  quite  substantially,  particularly  near  the  leading 

edge.  The  C  's  used  are  based  on  the  linear  formula  Reference  54)  and  a 
“o 

modified  linear-fit  of  Spreiter  s  nonlinear  formula  (Equation  (15), 

Reference  49);  again,  all  the  CD  coefficients  can  be  found  in  Appendix  C. 

In  Figures  21  and  22,  in-phase  and  out-of-phase  pressure  distributions  for 
the  same  given  conditions  (as  Figure  20)  are  compared  with  those  obtained 
using  Chadwick-Platzer ' s  MOC  method.  Observing  from  the  results  in 
Figures  19,  21,  and  22,  we  found  that,  unlike  the  transonic/subsonic 
results,  our  TLP2D  code  produces  generally  lower  in  pressure  magnitudes 
than  those  of  Chadwick-Platzer ' s.  However,  in  their  computation  procedure, 
the  steady  mean-flow  results  and  the  oscillatory  flow  results  are  all 
self-generated  and  the  former  result  is  truly  nonlinear  (based  on  Teipel, 
Reference  55).  It  would  be  interesting,  therefore,  to  check  against  these 
results  by  using  Teipel 's  mean-flow  input  in  the  TLP2D  program. 

5.4  Sonic  Results 

In  Section  3,  it  was  shown  that  the  sonic  kernel  function  (Equation  (3.25)) 
is  a  limiting  case  of  the  transonic  subsonic  kernel  function  (Equation  (3.15)) 
and  of  the  transonic/supersonic  kernel  function  (Equation  (3.22)).  Here, 
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a  numerical  study  was  performed  to  check  out  the  sonic  limit  for  these 
kernel  functions.  Figure  23  shows  that  although  the  kernel  functions  are 
continuous  throughout  the  sonic  limit,  the  LPP  method  (N  =  23  for  all  cases) 
produces  oscillatory  pressure  distributions  as  M®  approaches  one  (the 
parmeter  O  approaches  infinity).  Such  outcomes  were  found  by  Jordan  (Refer¬ 
ence  44),  using  the  pressure-mode  method,  and  was  subsequently  discussed, 
in  terms  of  the  wave  mechanism,  by  Landahl  (Reference  21).  This  then  led 
Jordan  to  be  in  favor  of  the  sonic  kernel  function  when  the  freestream  speed 
is  in  the  neighborhood  of  the  sonic  regime.  However,  the  objective  here  is 
a  rather  different  one.  We  attempt  to  establish  a  proper  mixed  kernel  func¬ 
tion  procedure  by  starting  out  with  a  case  of  sonic  flow  without  shock  wave 
(Figure  31),  Nevertheless,  in  so  doing,  one  must  check  out  the  sonic  kernel 
function  procedure  together  with  the  LPP  method.  This  then  constitutes 
the  following  calculation  results  using  the  sonic-flow  kernel  functions 
(Figures  24  through  30). 

Thus,  the  linear  sonic  flow  calculations  are  first  performed  for  a 
pitching  thin  airfoil  (flat  plate)  in  order  to  check  out  the  application  of 
the  TDL  code  at  M»  =  1.0.  Calculated  pressure  magnitudes  and  phase  angles 
are  compared  with  those  obtained  by  Stahara  and  Spreiter  (Reference  2b)  in 
Figures  24,  25,  and  26.  Excellent  agreement  can  be  seen  for  all  cases  con¬ 
sidered.  Figures  27  and  28  present  our  calculated  results  for  a  pitching 
Cuderley  airfoil  and  a  pitching  parabolic-arc  airfoil;  the  steady  mean-flow 
C  for  both  cases  are  taken  from  Reference  26  and  are  curve- fitted  in 
Appendix  C.  It  is  seen  that  the  present  results  in  pressure  magnitudes  are 
generally  lower  than  those  by  Reference  26  (similar  to  the  finding  of  our 
supersonic  pressure  magnitude  to  the  MOC  method).  Also,  considerable  dif¬ 
ferences  in  phase  angle  exist  at  various  reduced  frequencies  between  both 
cases.  We  believe  that  these  differences  could  be  caused  by  two  factors  in  the 
method  of  approach.  First,  the  unsteady  local  linearization  method  by  Stahara 
and  Spreiter  really  did  not  include  the  edge  condition.  Since  their  formula¬ 
tion  is  based  on  the  velocity  potential,  to  cope  with  such  a  condition  would 

\ 

require  the  Schwar tzchi Id  iterative  procedure.  Second,  in  their  analysis 
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(Reference  26),  the  transonic  parameters  appeared  in  the  kernel  potential  are 
all  expressed  in  terms  of  'receiving  point'  x  rather  than  the  'source-point' 

§.  This  amounts  to  ignoring  the  higher-order  flow  influences  due  to  the 
panels  upstream.  Figure  29  shows  the  comparison  of  the  results  between 
using  the  full  equation  (2.11),  and  using  the  approximate  equation,  (2.15). 

It  is  seen  that  the  deviation  in  pressure  distribution  is  slight  when  the 
reduced  frequency  is  low.  Figure  30  exhibits  a  comparison  of  all  previous 
sonic  methods  with  the  present  LPP  method  for  a  pitching  parabolic-arc  air¬ 
foil  (t  =  0.06)  at  a  constant  angle  of  attack.  Instead  of  using  the  steady 
mean-flow  Cp^  given  in  Appendix  C,  we  adopted  Dowell's  value  (aQ  =  0.45  and 
b0  =  0.585,  Reference  23),  which  amounts  to  X  =  0.585  and  T  =  0.45  for  our 
inputs.  Consequently,  the  lifting  pressure  plot  shows  that, once  again,  our 
value  is  of  the  lowest  magnitude,  but  it  follows  the  same  behavior  of  Goodman's 
theoretical  distribution. 

Finally,  in  Figure  31,  the  mixed  kernel  function  procedure  is  tested, 
using  a  nonlinear  sonic  flow  example.  For  simplicity,  the  Guderley  airfoil 
is  chosen,  in  which  the  sonic  point  is  located  exactly  at  40  percent  chord. 
Altogether,  33  panel  elements  are  employed,  with  15  equal-space  elements 
(5  percent  chord  each),  four  compact  elements  (occupying  5  percent  chord)  at 
the  leading  edge  and  14  compact  elements  (occupying  10  percent  chord)  in  the 
neighborhood  of  the  sonic  point.  In  fact,  the  panel  width  is  based  on  the 
criterion  that  each  panel  width  be  restricted  by  Ax  ^  tt/4  |c|  .  While  the 
transonic/subsonic  kernel  function  (Equation  (3.15))  is  applied  to  the  ele¬ 
ments  throughout  the  whole  airfoil,  the  transonic/supersonic  kernel  function 
(Equation  (3.22))  is  applied  only  to  those  downstream  of  the  sonic  point.  At 
the  sonic  point,  which  is  situated  between  subsonic  and  supersonic  panels,  the 
nonlinear  sonic  kernel  function  (Equation  (3.25))  is  used.  The  calculated 
result  is  not  altogether  satisfactory.  It  is  seen  that  MKF  results  follow 
closely  to  the  ones  by  the  sonic  kernel  function  method,  but  that  cusps  occur 
in  these  results  at  the  sonic  point.  At  this  stage  of  the  MKF  development, 
it  is  not  clear  whether  the  sonic  cusp  is  due  to  an  erroneous  kernel  func¬ 
tion  evaluation  or  due  to  an  inappropriate  paneling  scheme  in  the  neighborhood 
of  the  sonic  point.  Future  investigation  on  this  problem  is  mandatory  for  the 
development  of  the  MKF  method. 
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SECTION  6 


CONCLUSIONS 

A  transonic  kernel  function  method  for  unsteady  flow  calculations  has 
been  developed.  The  basic  formulation  of  the  transonic  acceleration  poten¬ 
tial  equation  was  derived  and  its  basic  solutions  in  two-  and  three-dimensional 
subsonic/ transonic/ super sonic  flow  regimes  were  obtained.  Hence,  the  two- 
dimensional  transonic/subsonic,  transonic/supersonic,  and  sonic  kernel  func¬ 
tions  were  then  established.  To  solve  for  the  downwash  integral  equation, 
two  discrete  element  methods  were  introduced;  these  are  the  Doublet  Lattice 
method  and  the  Linear  Pressure  Panel  method.  We  realized  that  the  former 
method  is  only  restricted  to  the  transonic/subsonic  (or  subcritical  and  purely 
subsonic)  flow  regime,  while  the  latter  method  has  the  unified  feature,  if  the 
control  point  is  properly  chosen,  for  subsonic/ transonic/supersonic  flow  cal¬ 
culations,  although  the  calculation  scheme  is  more  time-consuming  than  the 
former. 

To  demonstrate  the  transonic  LPP  method,  numerical  examples  were  pre¬ 
sented  for  simple  airfoil  geometries  in  pitching  motions  for  all  flow  regimes. 
For  thin  airfoil  cases  (i.e.,  for  purely  subsonic/sonic/super  sonic  cases), 
excellent  agreements  were  found  with  the  analytical  results  based  on  previous 
classical  theories.  For  transonic/subsonic  cases,  the  pressure  magnitude  is 
found  to  be  generally  larger  than  other  results,  while  for  transonic/supersonic 
and  nonlinear  sonic  cases,  it  is  found  to  be  lower  than  the  results  of  other 
theories.  Numerical  studies  revealed  that  discrepancies  and  sometimes  even 
numerical  divergence  in  the  present  transonic/subsonic  calculation  (e.g.,  Fig¬ 
ures  10  and  11)  are  largely  attributed  to  the  drastic  variation  of  the  tran¬ 
sonic  acceleration  parameters  and  F  in  the  neighborhood  of  the  leading  edge. 
Such  discrepancies,  in  principle,  can  be  greatly  improved,  however,  if  one 
adopts  Landahl's  general  formulation  of  phase  correction  method  (see  Equa¬ 
tions  (L.l)  and  (L.4)  of  the  Addendum).  Phase  angles  can  thus  be  corrected 
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accordingly,  however,  this  would  amount  to  developing  a  scheme  in  acoustic 
ray  tracing.  Lastly,  the  proposed  MKF  procedure  is  applied  to  a  nonlinear 
sonic  flow  case.  Close  agreements  are  seen  in  comparing  the  MKF  results 
with  the  nonlinear  sonic  kernel  function  results  except  at  the  sonic  point. 

It  is  not  clear  that  what  causes  the  undesirable  cusp  at  the  sonic  point  and 
how  to  smooth  itout  at  this  time. 

The  present  study  leads  to  the  following  concluding  remarks. 

1.  We  have  demonstrated  that  the  combination  of  the  transonic  kernel 
function  method  with  a  properly  chosen  discrete  element  method  such 
as  the  LPP  method  is  a  suitable  method  for  unified  treatment  of 
unsteady  transonic  flow  calculations. 

2.  Our  calculations  among  other  unsteady  flow  methods  assure  th»  dis¬ 
tinctive  departure  of  the  nonlinear  results  from  the  linear  ones  in 
all  transonic  flow  regimes  (not  to  mention  the  supercritical  regime). 

3.  The  present  study  indeed  paves  the  road  for  the  MKF  method  develop¬ 
ment  in  the  supercritical  regime.  It  also  opens  up  ways  for  the 
future  development  of  a  three-dimensional  LPP  method. 

In  view  of  the  encouraging  results  obtained  by  the  present  transonic  LPP 
method,  we  recommend  the  following  projects  for  further  development  and 
extension  of  the  present  work. 

1.  Improve  the  LPP  method  algorithm  so  that  the  TLP2D  code  can  be  made 
equally  cost-effective  as  the  TDL2D  code. 

2.  Generalize  the  TLP2D  code  to  include  three  degrees  of  freedom,  i.e., 
pitching,  plunging,  and  flapping  in  all  flow  regimes. 

3.  Incorporate  Landahl's  general  formulation  with  phase  corrections  into 
our  TLP2D  code. 

U.  Improve  the  MKF  procedure  for  smoothed  sonic  flow  cases  and  develop 
the  MKF  method  so  that  the  TLP2D  code  can  cope  with  supercritical 
flow  regime,  as  this  can  be  considered  as  the  most  important  step  in 
the  MKF  method.  In  this  regard,  two  steps  nu  suggested. 
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Computation  of  the  Eckhaus-Landahl  flow  model  (see  Appendix  B). 

Computation  of  the  oscillatory  supercritical  flow  with  embedded 
shock  waves. 


r 


The  basic  differential  equation  for  the  acceleration  potential,  t|r,  reads 
(corrected  and  valid  for  low  frequencies  only,  two-dimensional  case  considered) 


ft)  txx  +  tzz  +  A*x  +  Bit  =  0 


(L.l) 


where 


Po2  =  Poo2  -  h»2(y  +  iHx 


A  =  -  2ikMoo3  ( 1  +  _♦*)  -  2Moo3[(y  +  1)  -  fe®  <¥_;_l)]txx 


Poo2  =  1  -  Moo2 


For  Moo  close  to  unity  and  for  low  frequencies,  the  terms  underlined  may  be 
neglected.  Thus,  the  differences  between  this  result  and  the  earlier  one  by 
D.  Liu  are  unimportant. 


To  find  an  approximate  solution  valid  for  the  receding  wave  portion  we 


set 


♦  =  tie' 


Qf(x) 


(L.2) 


and  choose  df(x)  such  that  the  equation  for  \|r  has  only  0(1)  coefficients.  One 
finds  the  following  equation  for  ti : 


ft)2+lxx  +  tlzz  +  Aitix  +  Bit!  -  0 


(L.3) 


where 


Ai  =  A  +  2pb2</ 


Bi  -  B  +  Pq  a  +  Pb  («  )  +  Aar 

and  the  prime  superscript  is  defined  as  a  total  differentiation,  i.e., 
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We  set  f  =  fe2a.  Then,  with  ftj2  *  g^2  -  M2(y  +  1)  ,  we  have 

’  Ai  =  -21WW3  +  2(g0a)/  +  2f 

a 

B,  =  M»2[k3  -  <Y  +  D^-  21ky^xx]  +  -L.  [-2ikMco3  +  2(P0a)/] 


For  the  receding  wave,  we  anticipate  that  Of'  =  0(^j"3),  i.e.,  that  f  =  0(1). 
Hence,  in  order  for  the  coefficient  Bj  to  be  finite  in  the  limit  £^2  ■»  0  we 
must  have 


-2ikMoo3  +  (go2/  +  f  =  0 

or 

2ikMo>3 

a  =  -.MPe3)  +  f - dx 

Bo2 


With  this  substitution  we  may  set 


r  2ikM»2  hv 

*  -  fcfc-- 

where  is  one  solution  of  the  equation 


(L.4) 


@o2W  +  2ikH.a*ix  +  flzJ5  +  Moo2  [ka  -  aiky^Jh  -  0  (L.  5) 

The  solution  of  this  equation  may  be  expected  to  exhibit  slow  variations 
with  respect  to  x,  since  the  shortwave  component  should  be  described  by  the 
phase  factor  exp  {J'(2ikMoo2/@b2 )dx] .  Therefore,  the  first  term  may  be 
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neglected.  Also,  in  the  last  term,  2iky$Xx  could  be  neglected  compared  to 
k2  (for  low  frequencies,  both  k2  and  2iky4xx  may  be  neglected).  Note  the 
plus  sign  in  front  of  the  tylx-term.  Compare  this  with  Landahl's  result  in 
Symposium  Transsonicuum  I  (Reference  56).  The  result  for  i|f  thus  obtained  is 
only  valid  for  the  receding-wave  portion.  For  the  advancing-wave  portion, 
there  should  be  no  factor  0Q-2. 

We  shall  now  consider  a  recipe  for  constructing  a  uniformly  valid  approx¬ 
imation  for  the  pressure  kernel  starting  from  the  local  linearization  result, 
i.e.,  the  solution  obtained  by  assuming  a  uniform  flow  of  Mach  number  equal 
to  the  value  at  the  sonic  point.  Denote  this  solution  by  superscript  L.  The 
kernel  function  is  obtained  from  the  doublet  solution  and  from  the  above, 
one  can  prepare  the  following  approximation: 


*d 


CHQ(V)  (L)  i(0  -  0(L)) 
cn(v)  *d 


where  is  the  group  velocity  of  acoustic  waves  having  wave  front  at  an 
angle,  v,  to  the  x-axis  (see  Figure  LI),  c  ,  the  value  at  the  source  point, 
0,  the  phase  angle  of  the  acoustic  waves  along  a  particular  ray  connecting 
source  and  field  points,  and  0^L\  the  phase  of  the  local  solution. 

This  approximate  is  consistent  with  the  acoustic  ray  theory  (see 
Landahl's  article  in  Symposium  Transsonicuum  II,  Reference  53). 


Figure  LI.  Group  Velo¬ 
city  of 
Acoustic 
Waves 


c  =  ?*n=a-(l+4>  )sinv  +  4>  cosv 
n  x  z 

where  a  is  the  speed  of  sound. 
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For  the  two-dimensional  case,  only  the  angles  v  =  rr/2  and  -tt/2  appear,  so 
for  the  receding  wave 


cn  "  a  "  (1  +  V  ~  M»‘X  -  1  *  M»3(V  +  D«x 
(Notes  0Q2  =  1  -  M2  =  1  -  (l  +  ^x)2/a2  -  -  Mw2(Y  +  1)^) 


hence 


a  -  (1  +*x)  =  a2(a  +  1  +  *x)"1[eoo2  -  (y  +  1)^] 

“  ice.3  -  ew  +  ihxi  =  ieo3 

The  phase  angle  may  easily  be  calculated  in  the  two-dimensional  case.  During 
the  time  'dt'  the  acoustic  wave  propagates  a  distance  dx  =  c^dt  and  hence 


xo 


J* 


-JL- 


XO 


X 


2d£ 

0o3 


is  the  phase  shift  between  the  source  point,  Xq ,  and  a  field  point,  x. 

In  the  three-dimensional  case,  the  calculation  of  the  phase  shift 
becomes  more  complicated.  It  will  then  be  necessary  to  trace  the  acoustic 
rays  emanating  from  each  source  point,  and  the  wave  front  angles  along  the 
way.  In  the  ordinary  linearized  case,  the  rays  are  straight  lines,  but  in 
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the  transonic  case  they  would  be  curved,  particularly  for  the  receding  wave 
portions.  I  suggest  that  a  ray  tracing  program  be  worked  out  so  that  one 
can  see  how  rays  will  behave  for  a  particular  three-dimensional  wing  case. 
Possibly,  the  ray  curvature  may  be  small  so  that  one  could  approximate  the 
ray  path  by  a  straight  line  between  the  source  and  the  field  point,  in  which 
case  the  determination  of  the  phase  angle  simplifies  considerably. 


Going  back  to  the  two-dimensional  case,  the  kernel  is  obtained  from  the 
integral 

K  •  Urn  fJX  e-1W*  -  (x.jdxj 

Z-K)  -<*>  2 


lim  fj 

Z-*0  -co 


X  e-ik(x  -  Xl)  Cno  ^  <L>(Xl)eif0(xl>  *  e(L><xin 
cnvxi)  dz 


dxx  ] 


For  small  k, 
from  the  term 
celled  by  the 


which  is  the  main  concern  here,  the  main  phase  variation  comes 
exp(iQ),  the  phase  variation  in  the  linear  solution  being  can- 
term  exp(-i0^L^).  Now,  from  the  receding-wave  solution 


H  - 


Jc_ 

c 

n 


Hence,  integration  by  parts  of  the  above  approximate  kernel  gives 


K 


r  nO 

£  *<> 


(L)  i (0  -  e 
e 

z 


(L) 


) 


m 

J 


e-ik(x  -  xx)^ 


(L)(xl)ei^6(Xl) 


e(L)(x1)] 


dxj  } 
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Of  the  two  terms,  the  first  one  would  be  the  most  Important  one  for  low 
values  of  k.  If  we  now  compare  with  the  kernel  obtained  for  the  purely 
linear  case 


K  -  lim  *  <L)  c  fX  <>'ik(x  "  xi>*  1 

K  lir  *dz  •  %  J-e  *d,  dx^} 


z-*0 


we  find  that  a  good  approximation  should  be  obtained  by  taking  the  local 
linearized  kernel  function  for  the  local  Mach  number  at  the  source  point 
and  correcting  it  for  the  phase  angle.  Thus, 


Ktrans 


K(L)ei(0  -  0(L)) 


where  0  and  0^L^  refer  to  phase  angles  between  source  and  field  point  loca¬ 
tions.  Hence,  the  correction  involves  only  the  computation  of  phase  shifts 
using  ray  theory. 


In  order  to  carry  out  the  integrations  for  (denoted  by  AD^)  such  as 
Equations  (4.19),  (4.20),  and  (4.21),  it  is  necessary  to  first  curve- fit  the 
analytical  function  Fm(Xj  -  |)(m  =  0,1,2),  as  shown  in  Equations  (4.22)  and 
(4.23)  for  the  subsonic  case.  Notice  that  all  the  coefficients  used  (i.e., 

A,  B,  and  C)  in  F.^  are  not  those  defined  in  Equations  (2.11),  (2.15),  and  (2.16) 
in  the  text;  also,  6^  in  the  subsequent  formulas  should  be  distinguished  from 
Dirac's  delta  function,  6. 


As  the  subsonic  flow  expression  was  given  in  Equations  (4.22)  and  (4.23), 
we  shall  only  present  the  sonic-flow  formula  and  the  supersonic- flow  formula 
for  AD^. 


•  Sonic  Elements 


Fo(xj  -  5)  -  /Xj  -  5  KQ(Xj  -  £)  ,  x.  2  5 


Fitting  F0  as  a  parabola,  we  write 


Fq (Xj  -  l)  2*  A(x^  -  l)B  +  B(Xj  -  5)  +  C 


ADji  =  YogTO  j-  15A(a7/2  -  x.7/2)  +  2l(aA  -  B)(a6/2  -  x.6/2  ) 


where 


+  35(aB  -  C)(a3/2  -  Xj3/2)  +  105aC(a1/2  -  x^1/2)J 


a  =  Xj  -  S±-  «t 


(  A.  1 ) 


*i  " 


-«i  -  5i-x) 


for 


Si-i  s  5  s  5i»  x,  ^  5i 


Si  s  5  s  5i+l  »  x  ^  5  i+i 


5i+i  -  Si 


p 


Sonic  Element 


For  Xj  =  l/2(§i-i+  5i)  =  Xi.x, 

=  v/<5i  -  Si-l>/2  [Fo(xi-i-  §i_i )  +  4Fo(0)] 

For  Xj  =  (§i  +  |i+x)/2  =  xlf 


n/(?1-K  -  5i>/2  [2F0(xt  -  li)  +  3Fo(0)] 


Figure  A2.  Sonic  Elements  Showing 
Receiving  Points 
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•  Supersonic  Elements 


Fl(xj  -  5)  "  K^xj  -  5)  -  fl^xj  "5)  xj  *  l 

here,  6  is  Dirac's  delta  function.  For  *  5^* 


(Si  - 


[FX<xj  -  Si >  +  2F1(xj  -  —  Y1"1)]  (A*4) 


For  xj  *  Si-n, 


„  (§i+i  ”  §i7 

V  = — n?r — 

p 


[fx<xj  -  §i)  +  2Fj  ^xj 


5i  +  §i+i 


xj 51  ?i  Pi 


5i-i  xj  ?i  5 


(A.  5) 


X  2  l 

j  *1  +  1 


«i  «i+l  5 

Figure  A3.  Supersonic 
Element 

For  x j  =  (§i_i  +  Si ) / 2  =  xi„x, 


si  j  *i+l 


Figure  AA.  Supersonic  Element 

Showing  Receiving  Points 


(Si  -  Ci.i) 


[Fl(xi-l  *  ?i-x>  +  2F1(0)1 


xj  =  (5i  +  §i+i>/2  =  xi, 


<5i+i  -  to 


[5Fl  (Xi  -  10  +  AFj  (0)] 


(A. 6) 


(A.  7) 
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APPENDIX  B 

KERNEL  FUNCTION  FORMULATION  ACCORDING  TO  ECKHAUS-LANDAHL 

SHOCK- JUMP  MODEL 


AD-A080  227  NORTHROP  CORP  HAWTHORNE  CA  AIRCRAFT  GROUP  F/G  20/* 

TRANSONIC  KERNEL  FUNCTION  METHOD  FOR  UNSTEADY  FLOW  CALCULATIONS— ETC (U) 
MAY  79  D  D  LIU»  W  S  PI»  M  T  LANDAHL  F33615-78-C-3202 

UNCLASSIFIED  NOR-79-52  AFFDL-TR-79-3085  NL 


20=2 


AO 

A080227 


END 
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The  oscillatory  shock- jump  relation  based  on  Eckhaus-Landahl ' s 
one-dimensional  model  (Reference  21)  can  be  written,  in  terms  of  the  velocity 
potential,  <p,  as 


where 


Notice  that  the  LHS  of  Equation  (B.l)  represents  the  downstream  condition 
immediately  behind  the  shock  wave,  whereas  the  RHS  represents  the  upstream 
condition  immediately  ahead  of  it.  For  simplicity,  we  let  the  local  Mach 
numbers  in  Regions  1  and  2  be  a  uniform  supersonic  Mach  number  and  a  uniform 
subsonic  Mach  number.  The  mean  shock  jump  is,  of  course,  governed  by  the 
Prantdl  shock  relation.  Let  the  shock  point  be  xs,  and  hence  the  control 
points  in  the  panel  and  '+'  (see  Figure  Bl)  are  denoted  as  xg  and  xg  . 
Hence,  M  =  M,  and  M(x  +)  =  Mg  . 


Mean  shock 


Region  2 
<  1 


Region  1 
M,  >  1 


Aileron 


Figure  Bl.  Eckhaus-Landahl  Model  Showing  Notations 


91 


-r-v-rr-.-. 


According  to  the  usual  linearized  definition  of  acceleration  potential 
(i.e.,  dropping  the  <J>x  terms  in  Equation  (2.8)),  Equation  (B.l)  can  be 
recast  in  the  following  form,  after  a  di fferentation  in  the  z-direction,  i.e., 

tz+  -  i(k  +  Of)«Pz+  =  -*z"  +  i(k  -  a)cpz'  (B.2) 

Now,  if  we  define  an  operator  L0  as 

LoC-3"  «’kX  •  liml-  [•]  (B. 3) 

z-»0  dz 

then  L0  can  be  applied  to  Equation  (B.2)  and  it  yields 

K+  h  K+(Xs+,0)  =  {(a  -  k)K"(Xs",0)  -  i(L0[*z+]+  L0[*z'])}  (B.4) 


where 


X 

s 


X 

s 


-  5  ,  C  i  Xs 


+ 

X 

s 


5 


1  >  c  >  xs 


K'(xs',o) 


-ikX 

e 


The  last  two  terms  in  Equation  (B.4)  can  be  derived  from  Equations  (3.8)  and 
(3.9),  i.e., 

L°[*z*]s  x  _  x  + 

-  -cakMBe‘i<K,>,B"k)Xs+  •  H1(3,(xtXs+)/xs+  (B.5) 
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(B.6) 


-i(»c1Hl+k)X  ‘  r  .  _  -I 

-  Cle  8  •  [H(Xs")kM1Jl(KxX#  )/Xs  -  «na6,(Xs  )J0UlXs  )J 


where 


2n 

1 

ma  -  Mi"  - 

m 

in 

T 

{ 

8"  -  1  -  M 

kMj  /m* 

kMft/6* 

Consequently,  when  the  condition  at  x  ■  x  Is  obtained,  the  integrated  poten- 

s 

tlal  In  Region  2  can  be  written  as 


<Pq<X,z)  - 


-ikX 

e 


X  ikxo  ikX  + 

J  e  tdxo  +  e  <p(X  ,z) 

+ 

X 

s 


(B.7) 


Next,  let  us  introduce  another  operator,  L 0 [  •  1  s  lim  r-j  [.},  and  apply  this 

z-*0 

operator  to  Equation  (B.7).  After  some  rearrangement,  we  obtain  the  kernel 
function  formula  in  Region  2,  i.e., 


Kq(X,0>  -  K,(X,0)-  [k,(Xs+,0)  -  K+j 


-ik(x-xs+) 


(B.8) 


Notice  that  Kt  in  Equation  (B,8)  is  the  usual  subsonic  kernel  function  in 
the  absence  of  the  shock  wave. 
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For  the  application  of  the  above  kernel  function  formulation,  we  select 
the  aileron-buzz  model  proposed  by  Eckhaus  (Reference  52),  see  Figure  Bl. 

Thus,  the  shock  wave  is  placed  ahead  of  the  hinge  point.  As  we  have  assumed 
that  the  mean  shock  wave  is  one-dimensional,  which  is  infinite  in  extent, 
the  acoustic  waves  generated  by  the  oscillating  aileron  cannot  diffract 
toward  Region  1  upstream.  (Such  a  wave  diffraction  should  be  accounted  for 
in  the  two-dimensional  case,  as  it  was  shown  by  Tijdeman,  Reference  14. )  Ail 
the  upstream  influences  on  Region  2  from  the  supersonic  side  is  taken  care  of 
by  the  shock- jump  condition.  Hence,  in  terms  of  the  present  panel-method, 
the  above  consideration  amounts  to  no  further  influence  of  the  panels  betwen 
Region  l  and  Region  2  other  than  the  shock-jump  kernel  formula  relating  the 
adjacent  panels  and  '+'. 

For  a  calculation  example,  we  plan  to  compute  the  stability  boundary 
for  the  hinge  moment  (see  Figure  3,  Reference  52).  Since  our  formulation  is 
not  restricted  to  low-frequency  computation,  our  result  may  then  be  compared 
with  Lambourne's  experimental  data  (k  below  0.1). 

Although  Eckhaus'  formulation  is  fairly  general  in  the  frequency  range, 
his  calculation  scheme  is  subject  to  an  approximation,  which  restricts  the 
frequency  to  a  higher  range  than  0.1.  Furthermore,  as  the  present  MKF 
formulation  is  completely  flexible  in  the  steady  flow  input,  little  problem 
is  expected  to  generalize  the  aileron-buzz  model  to  one  with  a  nonuniform 
mean  flow.  Clearly,  studying  this  problem  will  certainly  better  our  under¬ 
standing  of  the  mixed  kernel  function  procedure  leading  to  a  general  scheme 
for  the  supercritical  flow  computations. 

Finally,  it  should  be  pointed  out  that  the  present  shock-jump  kernel 
function  formulation  only  amounts  to  one  which  satisfies  the  shock  jump 
condition,  Eq.  (B.l),  at  the  foot  of  the  shock,  i.e.,  at  0.  As  the 
assumed  moan  shock  wave  is  infinite  in  extent,  Eq.  (B.l)  should  be  satisfied 
formally  along  the  shock  in  z-direction.  Landahl  (Reference  57)  proposed  to  use 
Fourier  Transform  in  obtaining  a  general  two-dimensional  condition  for  a 
pulsating  source  and  has  obtained  solutions  in  the  upstream  and  the  down¬ 
stream  sides  of  the  shock  wave  in  the  transformed  plane.  Further  investi¬ 
gations  along  this  line  are  needed. 


0.05657  -0.11261 
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